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Multi-dimensional density of states provides a useful description of complex frustrated 
systems. Recent advances in Monte Carlo methods enable efficient calculation of the density 
of states and related quantities, which renew the interest in them. Here we calculate density 
of states on the plane (energy, magnetization) for an Ising Model with three-spin interactions 
on a random sparse network, which is a system of current interest both in physics of glassy 
systems and in the theory of error-correcting codes. Multicanonical Monte Carlo algorithm is 
successfully applied, and the shape of densities and its dependence on the degree of frustration 
is revealed. Efficiency of multicanonical Monte Carlo is also discussed with the shape of a 
projection of the distribution simulated by the algorithm. 



§1. Introduction 

Given a system of the energy E(x), the density of states of physical quantities 
Aq(x), Ai(x), ... is defined as 

D(A = a Q , Ax = a x , . . . ) = ^ j<5(A (a;) - a ) x 5(A x (x) - a x ) x • • • j 

where 5 is defined as usual: 5(a) = 1 if a = 0, else 5(a) = 0. Usually, one of the A^s' is 
the energy E(x) itself. When we set Aq(x) = E(x) and ignore other Ajs', it gives the 
original definition of a univariate "density of state" D(E), which is the normalization 
constant of micro canonical distribution P mc (x) = 5(E(x) — Eq)/D(Eq) 4 Hereafter 
we are mostly interested in bivariate cases where Aq(x) = E(x), A x (x) = M(x), 
where M(x) is an order parameter of a system, for example, magnetization of an 
Ising ferromagnet. 

Recently, an introduction of extended Monte Carlo methods, 1 ) especially multi- 
canonical Monte Carlo algorithm, makes calculation of multi-dimensional densi- 
ties of states a realistic choice for the analysis of complex probability distributions. 
There is a special interest in studying the cases with a first-order transition with 
latent heat with this approach. Existence of the latent heat indicates a rapid drop of 
the entropy of the system in ordering and suggests the difficulty in searching ground 
states. It will be interesting to investigate the nature of ordering process through the 
calculation of the multi-dimensional density of states. In this paper we study Ising 
models with three-spin interactions on random sparse networks ("Sourlas code") 4 ) 
and calculate density of the states and related quantities by multicanonical Monte 
Carlo. 
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§2. Multicanonical Monte Carlo 

Details of multicanonical Monte Carlo algorithm are found in the references. 1 ^' 2 ) 
The essence of the algorithm is dynamical Monte Carlo sampling with a weight 
1/D(E), where D(E) is an approximation of D{E). It causes a nearly uniform 
distribution of the energy E within an interval. The density D(E) is estimated by 
repeated "preliminary runs" of the simulation (a simple method called "entropic sam- 
pling" 5 ) is used here). The bivariate density D(E,M) is obtained by log D(E, M) = 
log W(E, M) + \og(D(E)) + const. j where W(E, M) is the frequencies of the values 
(E,M) appeared in a simulation with the weight \/D(E). We can also reconstruct 
a conditional density D{M\E) by log D(M\E) = log D(E,M) - log D(E, M). 
D(E\M) is proportional to the projection on (E,M) plane of the distribution sam- 
pled by multicanonical Monte Carlo with an ideal weight D(E) = D(E). 

The advantage of multicanonical Monte Carlo is that relaxation in small E (and 
also very high E) region is greatly facilitated compared with canonical or micro- 
canonical simulations. An interesting question is when multicanonical Monte Carlo 
performs better than the other methods with extended ensembles, parallel temper- 
ing 1 ) and simulated tempering. 1 '> 6 ) It is usually believed that multicanonical Monte 
Carlo performs better in cases with first-order transitions with latent heat, because 
multicanonical ensemble can contain states that rarely appear in a canonical dis- 
tribution of any temperature, and relaxation is speed up with a vanishing critical 
nucleus size. This picture, however, seems not to be tested enough. We will discuss 
it in the later section with our example. 

§3. Ising Model with three-spin interactions 

Let us consider an Ising model with three-spin interactions: 

Pp(x) oc exp \/3 ^2 J ijk x 

\ (ijk)eG J 

Here Yl(ijk)eG denotes the summation on the edges of a graph G. The energy and the 
magnetization is defined as E = Y^Ujk)eG JijkXiXjXk, and M = Yli x ii respectively. 
Several different definitions of a random network G is possible. Here we consider a 
simplest case to simulate: The probability of the presence of the edge is independent 
and takes the same value for all pairs of the vertices on G. Total number of edges 
are fixed to a given number I, but the number K of the edges that contains a vertex 
is not fixed. Hereafter we set Jijk = ±1. 

An example of the results are shown in Figs. ^ and |2] for a ferromagnetic model 
(Jijk = 1) °f the si ze N = 100 on a random graph G with I = 200 edges. Density 
of states D(E,M) and the projection D(E\M) of multicanonical density are plotted 
in Figs. ^ and |21 The vertical and horizontal axes corresponds to the energy E and 
the magnetization M. In Figs. Q and |2J an interval of the contours corresponds to 
the change of the density by the factors 10 4 and 10 ' 5 ~ 3.2, respectively. 
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Fig. 1. Density of States: log D(E,M) Fig. 2. Multicanonical Density: log D(M\E) 

In Fig. ^ states around (E,M) = ±(200, 100) at the ends of "arm"s observed in 
both sides of the body correspond to the ordered states, while the peak at the center 
corresponds to the paramagnetic state. Point symmetry of the density corresponds 
to Si — > — Si, (E, M) — > (— E, —M) symmetry of the model. In Fig.|2 three peaks of 
the density are observed, one of which is the paramagnetic state at the center, and 
the others correspond to ordered states. Multiple peaks indicate that metastabilty 
appears even with multicanonical Monte Carlo. 
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Fig. 3. Distribution of (E, M ) at = 0.4. Fig. 4. Distribution of E at (3 = 0.4. 

In Figs. [Handel projections of a canonical distribution Pp{x) ((3 = 0.4) on the 
(E, M) plane and on the E plane are shown, respectively. In Fig. |3J the contour 
interval and axes are the same as those in Fig. [21 Two peaks are again seen in both 
figures, which agree with an existence of the first-order transition with latent heat. 4 ) 

The changes of the shape of density of states when we add frustration are shown 
in Figs. 15171 Models with N = 100 and I = 500 are investigated for the ratio p 
of negative bonds 0, 0.1, and 0.36. The contour interval and axes are the same as 
those in Fig. ^ Increasing the concentration of frustration, the maximum value of 
magnetization that appears with a probability higher than a threshold reduces and 
the "arm"s shrink. At the value p = 0.36, there are no more evident arms, which 
indicate ordered states do not exist at any temperature. 
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Fig . 5 . p = Fig. 6. p = 0.1 Fig. 7. P = 0.36 

For each of these examples, 1.0 ~ 1.5 x 10 9 MCS are required for the estimation 
of D(E) and ~ 10 s MCS are used for the computation of the physical quantities. 
Typical CPU times are 10 ~ 15 hours with a 3.2GHz Pentium4 chip. 



§4. Discussion and Future Problems 



The result shown in Fig. give rise to the possibility that (univariate) multi- 
canonical Monte Carlo may not have a clear advantage over parallel tempering in 
this example. A natural question is that how multicanonical density looks like for 
other systems with first-order transition with latent heat. Our preliminary results on 
10-states Potts model on the square lattice indicate that multicanonical density pro- 
jected on a (energy, order parameter) plane shows smoother behavior compare with 
those of Ising models with three-spin interactions, although a saddle point is still 
observed. Detailed comparison is, however, left for future studies, as well as com- 
putation of a thermodynamic limit of the densities and direct comparison between 
multicanonical and parallel tempering algorithm. 

Another interesting issue is the use of two-dimensional multicanonical Monte 
Carlo 3 ) with the weight 1/D(E, M) that approximates inverse of the bivariate den- 
sity D(E, M). It enables efficient sampling of small D(E, M) regions, and, hopefully, 
further enhance the mixing. Our experience, however, shows that it requires con- 
siderably many iterations for the estimate of the bivariate D(E,M) by entropic 
sampling. Some improvements in the learning stage should be introduced in future 
studies. 
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